Wannier-function approach to spin excitations in solids 
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We present a computational scheme to study spin excitations in magnetic materials from first 
principles. The central quantity is the transverse spin susceptibility, from which the complete 
excitation spectrum, including single-particle spin-flip Stoner excitations and collective spin-wave 
modes, can be obtained. The susceptibility is derived from many-body perturbation theory and 
includes dynamic correlation through a summation over ladder diagrams that describe the coupling 
of electrons and holes with opposite spins. In contrast to earlier studies, we do not use a model 
potential with adjustable parameters for the electron-hole interaction but employ the random-phase 
approximation. To reduce the numerical cost for the calculation of the four-point scattering matrix 
we perform a projection onto maximally localized Wannier functions, which allows us to truncate 
the matrix efficiently by exploiting the short spatial range of electronic correlation in the partially 
filled d or f orbitals. Our implementation is based on the full-potential linearized augmented-plane- 
wave (FLAPW) method. Starting from a ground-state calculation within the local-spin-density 
approximation (LSD A), we first analyze the matrix elements of the screened Coulomb potential 
in the Wannier basis for the 3d transition-metal series. In particular, we discuss the differences 
between a constrained nonmagnetic and a proper spin-polarized treatment for the ferromagnets Fe, 
Co, and Ni. The spectrum of single-particle and collective spin excitations in fee Ni is then studied 
in detail. The calculated spin-wave dispersion is in good overall agreement with experimental data 
and contains both an acoustic and an optical branch for intermediate wave vectors along the [1 0] 
direction. In addition, we find evidence for a similar double-peak structure in the spectral function 
along the [1 1 1] direction. To investigate the influence of static correlation we finally consider 
LSDA+t/ as an alternative starting point and show that, together with an improved description of 
the Fermi surface, it yields a more accurate quantitative value for the spin-wave stiffness constant, 
which is overestimated in the LSDA. 

PACS numbers: 71.15.Qe, 71.45.Gm, 75.30.Ds, 71.20.Be 



I. INTRODUCTION 



Spin excitations in solids are of fundamental interest 
for a wide variety of phenomena. For example, in mag- 
netic materials at low temperatures collective spin excita- 
tions, so-called spin waves or magnons, 1 leave a mark on 
the transport, dynamical, and thermodynamical proper- 
ties. The spin waves contribute to the specific heat with 
a T 3 / 2 term in addition to the T 3 term from phonon 
excitations. The low-temperature behavior of the mag- 
netization in three-dimensional magnetic solids is also 
dominated by spin- wave excitations: In ferromagnets the 
magnetization drops as T 3 / 2 , while the sublattice mag- 
netization in antiferromagnets obeys a T 2 law^ In low- 
dimensional systems spin-wave excitations even destroy 
the long-range magnetic order completely at any finite 
temperature T in the absence of a magnetic anisotropy. 3 
As the temperature increases, additional single-particle 
spin-flip processes, the so-called Stoner excitations, take 
place, which further contribute to the temperature vari- 
ation of the magnetization and cause a damping of the 
spin-wave modes. Spin waves can also couple to charge 
excitations around the Fermi level. Such interactions lead 
to a pronounced energy renormalization in the quasipar- 
ticle band dispersion,- or control the spin-dependent in- 
elastic mean free path of hot electrons in ferromagnets 



Another interesting phenomenon is high-temperature su- 
perconductivity, in which spin waves have been pro- 
posed as a possible mediator for the attractive electron- 
electron interaction^^ This scenario was recently con- 
firmed by infrared spectroscopy and neutron scattering 
measurements for the high-temperature superconductor 
YBa 2 Cu 3 6 . 92 Ai^ 

Spin excitations in solids are also of central interest in 
the field of spintronics. The writing process of magnetic 
information in a giant magneto resistance or tunnel mag- 
neto resistance device is closely related to the rotation of 
the magnetization, a process generating and radiating 
spin waves at all wave lengths, whose damping rate is 
an important parameter determining the writing timeJ^ 
Spin waves are also generated during the reading process, 
as hot electrons impinge at the interfaces between the in- 
sulating barrier and the magnetic electrodes of magnetic 
tunnel junctions, which causes a reduction of the mag- 
neto resistance^ 

The properties and physics of spin waves evidently 
comprise an unusually rich area of research. A lot of 
information about the spin dynamics in solids can be ob- 
tained from the magnetic response function (or dynami- 
cal spin susceptibility). The spectrum of magnetic exci- 
tations corresponds to the poles of the response function 
and can be directly compared with experiments like in- 
elastic neutron scattering. In this way it provides insight 
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into the nature of the exchange coupling and the com- 
plex magnetic order. The magnetic response function is 
thus a central quantity for the theoretical description of 
magnetic materials. 

So far most theoretical studies of magnetic excitations 
in solids were based on an adiabatic treatment of the 
spin degrees of freedom in which the slow motion of the 
magnetic moments and the fast motion of the electrons 
are separated by mapping the complex itinerant electron 
problem onto the classical Heisenberg Hamiltonian with 
exchange parameters obtained from constrained density- 
functional theory (DFT) calculations J£r— Within this 
approach the spin-wave excitations can be calculated ef- 
ficiently, whereas the single-particle Stoner excitations 
are neglected. Furthermore, the spin- wave life times are 
not accessible. From a fundamental point of view, the 
Heisenberg model is justified only for local-moment sys- 
tems like insulators and rare earths, which possess well- 
defined spin-wave modes over the entire Brillouin zone, 
so that the adiabatic approximation is reliable. For 
itinerant-electron magnets the adiabatic approximation 
yields reasonable results only in the long- wave-length re- 
gion, i.e., for small wave vectors. For short wave lengths 
the discrepancy with experiments can be large, however. 
For example, the multiple branches in the spin- wave dis- 
persion of 3d ferromagnets cannot be capturedi 23 ' 24 

First-principles calculations of the magnetic response 
function using realistic energy bands and wave func- 
tions are very rare. Initial attempts by Cooke el al*£~— 
within the framework of many-body perturbation theory 
(MBPT) were based on a tight-binding description of the 
electronic energy bands. These authors studied the spin- 
wave dispersion of 3 c? ferromagnets and obtained reason- 
able agreement with experiments; in particular, the op- 
tical branch in the spin-wave dispersion of fee Ni was 
correctly described.— Using a similar approach, Mills 
and co-workers 2 ^— carried out extensive calculations 
to explore the spin dynamics in ultrathin ferromagnetic 
films on nonmagnetic substrates. Recently more realistic 
treatments of the spin- wave spectra in 3 c? ferromagnets 
were reported by Savrasov^i and Buczek et al? 2 within 
time-dependent density-functional theory (TDDFT) and 
by Karlsson and Aryasetiawan within MBPT. 33 However, 
in the latter work the authors used a simplified model po- 
tential with an adjustable parameter to estimate the ma- 
trix elements of the screened Coulomb interaction. Un- 
der these circumstances both TDDFT and MBPT ap- 
pear to give similar results for the spin-wave dispersions 
of Fe and Ni. The results for Fe are in good agreement 
with available experimental data, while for Ni the optical 
branch is too high in energy, which was attributed to the 
overestimation of the exchange splitting of Ni within the 
underlying LSDA^ 3 - 

The aim of the present work is to develop a practi- 
cal computational scheme to study excitation spectra of 
magnetic materials from first principles. The magnetic 
response function is calculated within a many-body con- 
text following the formalism given in Ref.l34l. To study 



collective spin- wave excitations we include vertex correc- 
tions in the form of ladder diagrams, which describe the 
coupling of electrons and holes with opposite spins via 
the screened Coulomb interaction. In analogy to the 
T-matrix that describes the particle-particle scattering 
channel, here we use the same term for the electron-hole 
channel in agreement with the definition of Strinatij^ 5 . In 
contrast to earlier treatments, the matrix elements of the 
screened Coulomb potential are calculated entirely from 
first principles. In order to reduce the numerical cost for 
the calculation of the four-point T-matrix we exploit a 
transformation to maximally localized Wannier functions 
(MLWFs), which provide a more efficient basis to study 
local correlations than extended Bloch states^ - — This 
use of localized orbitals makes our scheme very efficient 
for complex magnetic materials with many atoms per 
unit cell. Our implementation is based on the FLAPW 
method. In the following we first calculate the matrix 
elements of the Coulomb potential for the 3 c? transition- 
metal series in the Wannier basis and perform extensive 
convergence tests. The magnetic excitations in fee Ni are 
then studied in detail based on the LSDA and LSDA+fJ 
methods. We find that both approaches yield qualita- 
tively similar results for the spin- wave spectra and overall 
dispersion. However, the static correlation effects seem to 
be important for the spin- wave stiffness constant, which 
is overestimated within LSDA. In contrast to some pre- 
vious theoretical studies ; 25 ' 33 our calculations clearly in- 
dicate the existence of an optical branch in the spin- wave 
dispersion curve of Ni along the [1 1 1] in addition to that 
along the [1 0] direction in the Brillouin zone. In gen- 
eral, the obtained results are in good agreement with 
available experimental data. 

This paper is organized as follows. In Sec. [II] we 
describe the computational method. Section llHl con- 
tains the results for the matrix elements of the screened 
Coulomb potential for the 3d transition metals. In 
Sec. lIVI we present results for the magnetic excitations 
in fee Ni together with a detailed discussion. In Sec.lVl 
we summarize our conclusions and give an outlook. Un- 
less otherwise indicated, Hartree atomic units are used 
throughout. 

II. COMPUTATIONAL METHOD 

A. Magnetic response function 

The time-ordered magnetic response function (or dy- 
namical spin susceptibility) is given in real space by the 
correlation function 

R«(l,2) = -i(T[a i (l),a j m, (1) 

where T is the time-ordering operator and <t*(1) are the 
spin-density operators with i S {x, y, z, — , +}, where 
— and + correspond to the spin annihilation (<r~ = 
a x — ia y ) and creation (<r + = a x + ia v ) operator, re- 
spectively. For simplicity we use the short-hand notation 
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1 = (ri,ti). The expectation value of with respect 
to the many-body ground state is given by 



(<7 i (l)> = -i5>^G a/3 (l,l+) 

a,0 



(2) 



with the Pauli spin matrices er 4 , the single-particle Green 
function G, and the spin indices a and j3. The notation 
1 + indicates that the time variable is increased by an 
infinitesimal to ensure the proper time ordering > t\. 
The magnetic response function can be obtained from 
the spin density by the functional derivative 



7^(1,2) 



(3) 



where i and j correspond to the components of the mag- 
netization and the magnetic-held vector, respectively. 
The latter incorporates a factor g/is/2, where the /j,b 
denotes the Bohr magneton and g the electron (/-factor, 
so that the Zeeman term in the Hamiltonian operator 
takes the form +B • a. 



B. T-matrix approximation 

The single-particle Green function in Eq. @ is given 
by the Dyson equation 

G af} (l,2) = G°(1,2)<W + E // G °( 1 ' 3 ) 



(4) 



xE a7 (3,4)G 7/3 (4,2)d3d4, 



where G° is the spin-diagonal Green function of the non- 
interacting Hartree system and S the nonlocal dynamic 
self-energy, which incorporates all exchange-correlation 
effects. With the identity 

SG a p(l,3) 
8B 3 {2) 

= "E J J G a7 (l,4)-^pG^(5,3)d4d5(5) 

one can rewrite the magnetic response function ([3]) as 
i?^'(l,2) 

k^(2 - 3) 



-i E a U J J G ai {l, 3) 
<5£ 7 a(3,4) 



x5(3-4) + 



^•(2) 



G^(4,l+)d3d4. (6) 



In this work we use the GW approximation for the self- 
E 7S (3,4) = *G 7ff (3,4)W(3 > 4). (7) 
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energy-^ 



The functional derivative of the self-energy with respect 
to the external magnetic held is then given by 



^(3,4) 
8B 3 {2) 



5W{Z,A) 
X ^B~W 



iG 7(5 (3,4) 



(8) 



For systems with a collinear magnetic ground state only 
the hrst term on the right-hand side yields a nonzero con- 
tribution to the magnetic response function^ Further- 
more, in this case the Green function is diagonal in spin 
space and can be written as G a/ g(l,2) = G a (l, 2)8 a p. 
The dynamically screened Coulomb potential W(3,4) is 
given by 

W(3,4) =u(3,4) + J J v(3,5)P(5, 6)1^(6, 4) d5d6, (9) 

where u(3,4) = 8(t% — £4)/^ — r±\ is the bare Coulomb 
potential and P(5, 6) the polarizability in the random- 
phase approximation (RPA). The latter is expressed by 



P(5,6) = -E^ QQ (5>6;6,5), 



(10) 



where the kernel K is defined as 

K<*P(1, 3; 4, 2) = iG a (l, 3)G^(4, 2+) . (11) 
After collecting all terms we obtain 

i?«(l,2) = -^^[^(1,2;2,1) 

a,f) 

+2/^(1,2; 2,1)] (12) 

for the magnetic response function. The second contri- 
bution is given by 

L af3 {l, 2; 2, 1) = JJJJ K a P(l, 3; 4, l)T a/3 (3, 5; 6, 4) 

X# a/3 (5,2;2,6)c23d4c25<26, (13) 

where the T-matrix obeys the Bethe-Salpeter equation 

T a/3 (1,3;4,2) 
= W(l,2)5(l-3)£(2-4) + W(l,2) 

x yy > X Q ' 3 (l,5;6,2)T Q ' 3 (5,3;4,6)d5d6. (14) 

The hrst term in Eq. (TT21 represents the response of the 
noninteracting system, i.e., the Kohn-Sham spin sus- 
ceptibility. The second term contains the T-matrix, 
which describes dynamic correlation in the form of re- 
peated scattering events of particle-hole pairs with oppo- 
site spins and is responsible for the occurrence of collec- 
tive spin-wave excitations. The Feynman diagrams for 
the magnetic response function R and the T-matrix are 
displayed in Fig.[TJ 
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FIG. 1: Diagrammatic representation of (a) the magnetic 
response function and (b) the T-matrix. 



For a practical evaluation of the magnetic response 
function we replace the full renormalized Green func- 
tion G Q (1,2) by an appropriate mean- field approxima- 
tion (LSDA, LSDA+f/, etc.). Moreover, we employ an 
instantaneous interaction of the form W(ri, r 2 ; t\ — t 2 ) ~ 
W(r 1 ,r 2 )5(* 1 -t 2 )withW(r 1 ,r 2 ) = W(r u r 2 ; r) dr 
in Eq. ([T4"]) . This is justified by the observation that the 
matrix elements of W(ri, r 2 ) do not vary strongly in the 
low- frequency region (w < 1 eV) , in which the spin- wave 
excitations occuri 55 ! 57 Under these circumstances it is 
sufficient to calculate the kernel ([TT]) for £ 2 = ti and 
ti = £3. In fact, due to time translation symmetry the 
resulting expression depends only on the difference t%— ts 
and is given by the Fourier transform of 



ir^(ri,r 3 ;r4,r 2 ; W ) 



occ unocc / a / \ a* t \ P f \ 0* f \ 



2. ^ v 

k , m k' ,rn' 



^ + ( £ k'm' 



"km 

/3* 



Vk'm'( r l)^k*m'( r 3)Vkm( r 4)¥'km( r 2) 



,) + »* 



(15) 



where <^£ m (r) and e£ m are the LSDA (LSDA+t/) eigen- 
states and eigenvalues, respectively. 



C. Implementation in the Wannier basis 

In order to reduce the numerical cost for the calcu- 
lation of the four-point kernel K we exploit a transfor- 
mation to maximally localized Wannier functions, which 
allows us to efficiently truncate the matrix in real space. 
The generalized Wannier functions w" R (r) with orbital 
index n and spin a at the site R arc defined as Fourier 
transforms of the Bloch states VkmO") according to 



»Sn(r) = ^E^^^E^-VkmW 

k m 

l£ e -*-R <n(r)j (16) 



where N is the number of discrete k points in the full Bril- 
louin zone and Umn^ denote the transformation matrices. 



The latter are determined by minimizing the spread 

= E («ok 2 Ko> - «o|r|< ) 2 ) , (17) 



where the sum runs over all Wannier functions. We em- 
ploy the algorithm for minimizing the spread initially 
proposed by Marzari and Vanderbilt2& for isolated groups 
of bands and later extended to entangled energy bands <21 
The matrix elements of the screened Coulomb poten- 
tial in the MLWF basis are given by 

W a/3 

niRinsRs ;n4R,4?i2R.2 

= J J <: Rl (r)< 3 R 3 ( r M r > r/ ) 

x<Vr')<R 2 (rVrdV. (18) 

The screened potential W(r, r') itself is calculated within 
the RPA using the mixed product basis4£~— As we use 
only the on-site matrix elements of W, we set Ri = R 2 = 
R 2 = R4 and eventually obtain 



a/3 



n 1 n 3 -,n4n 2 



N 3 E E E 

k,qi,q2 I,J rri! ,m 2 ,m 3 ,m4 



t ja (k+qi )* TT (k+q 2 ) 
u mini m-zTi^ 



X CCins * (Vk+di mi I Vq x m 3 Mfk) 

x(M JU |W(r,rO|M Jk )(M Jk ^ 2m J^ +q2m2 ) , (19) 

where Mjk (M/k) are the biorthogonal basis functions of 
the mixed product basis, which satisfy the relation 



7,k 



lAOkXMjkl = 1 



(20) 



in the Hilbert space of the wave-function products. 

The next step is the calculation of the kernel K . The 
projection of the Bloch states onto the Wannier orbitals 
yields 

/ <k(rM m (r) d 3 r = U^e** . (21) 

If we perform a lattice Fourier transformation, then 
Eq. (fTS"|) takes the form 



^nfn 3 ;n 4 n 2 (q, w ) 

1 occ unocc , TT a(\i)jja(Vi)* 1 -jfS(^+<i) T jP(^+ci)* 

1 \ - \ - \ - / Umni U m n 3 <->rn'm U m 'ni 



k m m' 



uj + (e 



— e a ) — iS 

k+qm' km/ 



7 -ra(k+q)* r .c t (k+q) r ./3(k)* 7 -./3(k) 
u m'ni u m>n 3 ° mn i u rnn 2 

W - Ok+qm< ~ 4m) + iS 



(22) 



in the Wannier basis. Instead of a direct evaluation of 
this expression, we first calculate the corresponding spec- 
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tral function 

1 occ unocc / 

_\ \^ \ ' f 7- / a(k+q)*rra(k+q) / -r/3(k)* 
Y / j / j / j \ u m'n 1 u m>n 3 mn 4 
k m m' ^ 

■xnP( k )x(t.) — f? -t- f 13 ) - r/ Q ( k )r/ Q ( k )* 

AL/ mri 2 tl, > u; fc k+qm' ^ fc krn/ ^mm ^ranj 

x«: q) ^ k „: q) ^(-+^ +qm .- e g m )) , (23) 

which equals the probability distribution for spin-flip 
transitions between occupied and unoccupied states with 
the energy and momentum difference uj and q. Once the 
spectral function is known, we use a Hilbert transforma- 
tion to calculate the kernel 



l>oo ca 



+ ^fn 3 ;„ 4 » 2 (q^)sgnM,(24) 

where P indicates the Cauchy principal value. 

With the kernel K and the screened Coulomb potential 
W we can construct the T-matrix according to Eq. (Tl4l) . 
which in the MLWF basis takes the form 



^nin3;n47i2 (q> ^) 



"tiin 3 ;n 4 ti2 1 / y 7li 

xT„f„ 3; „ 4 „ 8 (q, W ) 



(25) 



and can be solved by a matrix inversion for a set of q 
and w values. Finally, the magnetic response function is 
given by 

R«(q,u,) = -EE E 

a, /3 k 711,712,713,714 

x (q|<k+ q wf' k ) [^"f„ 3 ;„4n 2 (q> W ) 

«k +q <* k |q> (26) 

with 

<q|< lk+q <* k ) = / e- Jq - r < k+q (r)<* k (r)d 3 r. (27) 

The tilde denotes the orthonormalized products of Wan- 
nier functions. Although the Wannier functions them- 
selves form an orthonormal basis set with (w" k |u/^, k ,) = 
NSnn'Skk'Sap, their products do not satisfy this or- 
thonormality condition. Therefore, we explicitly or- 
thonormalize the products according to 

Kik+q*l*k) 

= E [o(q,k)- 1 / 2 ]:f n2i „ 3 „ 4 |< k+q <k) I (28) 

713,714 

where the overlap matrix is defined as 

0^ n2 ,n s nM> k ) = Wik+q^'kKk+q^k) ■ ( 2 9) 



In practice, this orthonormalization can be performed in 
the final step of the calculation, i.e., in the projection of 
the magnetic response function R onto plane waves. 

The spin-wave spectra are obtained from the imagi- 
nary part of the transverse magnetic response function 
R h (q, uj), which exhibits peaks at the spin- wave ener- 
gies corresponding to the wave vector q. The half-width 
of a peak is inversely proportional to the life time of the 
excitation. 



D. Computational details 

All ground-state calculations are carried out using the 
FLAPW method as implemented in the FLEUR code,— 
initially within the LSDA for the exchange-correlation 
potential^ We use 4.5 bohr -1 as a cutoff for the plane 
waves and l cu t = 10 for the angular momentum for all 
3d transition metals under consideration. In addition, 
the LSDA+C/ method with U = 1.9 eV and J = 1.2 cV 
is employed to reveal the correlation effects on the mag- 
netic excitation spectra of Ni. In practice, the U and J 
values can either be chosen as empirical parameters or 
obtained from from first-principles calculations by em- 
ploying methods like constrained LSDA 54 or constrained 
RPA, 5 — ~— in which the screening due to the 3d electrons 
is excluded. However, due to s-d hybridization in the 
3d transition metals the constrained RPA yields differ- 
ent values for U and J depending on the procedure used 
for excluding 3d-3d transitions in the polarization func- 
tion. For example, Miyake et al^ found 3.7 eV for U in 
fee Ni if the s-d hybridization is switched off, whereas U 
reduces to 2.8 eV if the s-d hybridization is retained. On 
the other hand, our calculations showed that using U and 
J values from the constrained LSDA or from constrained 
RPA within the LSDA+C/ scheme yields unsatisfactory 
results for the magnetic moment, exchange splitting and 
spin-wave dispersion of fee Ni compared to experiments. 
For this reason we use the empirical values for U and J 
given above, which improve the Fermi surface and do not 
change the magnetic moment substantially^ Further- 
more, these values yield the correct magnetic anisotropy 
energy and direction of the magnetization. 

The MLWFs are constructed with the Wannier90 
code^ which was recently interfaced to the FLAPW 
method*^ The screened Coulomb potential in the RPA 
is calculated with the SPEX cod o 50 ' 51 using the mixed 
product basis and then projected onto the MLWF basis. 
The total number of functions in the mixed product ba- 
sis is 180-200, and 95-100 unoccupied states are included 
in the calculation of the polarizability. Finally, we note 
that although the MLWFs provide a minimal basis set 
for the construction of the T-matrix, the computational 
time scales as the fourth power of the number of Wan- 
nier functions. The most expensive part in our scheme 
is the calculation of the kernel K, because it requires 
a large number of k points for proper convergence (see 
Sect.lIVp. In contrast, the screened Coulomb potential 
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W is less sensitive to the k-point sampling, i.e., it is al- 
ready converged for a substantially smaller number of k 
points. 

III. MATRIX ELEMENTS OF THE COULOMB 
POTENTIAL 

As a first step we calculate the matrix elements of 
the screened Coulomb potential W for the series of 3d 
transition metals, because these matrix elements are 
a crucial ingredient for the construction of the mag- 
netic response function. In previous treatments of spin 
waves in a many-body context the Coulomb interaction 
was chosen either as a simple Hubbard-type U param- 
eter or as a model potential with an adjustable range 
parameter . 24 ' 25 ' 33 Recent ab initio studies of the bare and 
the screened Coulomb interaction in 3d transition metals 
focused only on the nonmagnetic (NM) statej 55 ' 56 ' 61 " — 
Here we present a detailed study of the matrix elements 
in the MLWF basis for the proper ferromagnetic (FM) 
state of the 3d transition metals Fe, Co, and Ni. For 
comparison with previous works the NM states of these 
three elements and the rest of the 3d series are consid- 
ered, too. We focus especially on bcc Fe to investigate 
the effect of the exchange splitting on the Coulomb ma- 
trix elements, because bcc Fe has the largest exchange 
splitting among the 3d ferromagnets. Previous studies 
showed that, similar to insulators, in 3d transition met- 
als the Wannier functions with d character are exponen- 
tially localized^ 3 - The correlation effects hence take place 
predominantly within the same atomic site£l Our cal- 
culations confirm these findings. As the off-site matrix 
elements of the screened Coulomb potential are at least 
two orders of magnitude smaller than the on-site ones, we 
only consider the latter. The strong localization of the 
3d orbitals can be seen in Fig. [3J where we present e g -like 
(3d 3z 2_ r 2 and 3d x 2_ y 2) MLWFs for bcc Fe. The isosur- 
face corresponds to 10% of the maximum amplitude of 
the MLWFs. 

To begin with, we define the average on-site diagonal 
and off-diagonal matrix elements of the direct and ex- 
change Coulomb potential as 

(3d) 

U=zY, W ™mn, ( 3 °) 



(3d) 

V'=2Q E Cl;n„ (31) 
m, n 
(m ^ n) 

(3d) 

J =^ E W rnt,n m - (32) 
m, n 
(m ^ n) 




FIG. 2: (Color online) e g -like (3d 3z 2_ r 2 and 3d x 2_ y 2) maxi- 
mally localized Wannier functions for bcc Fe. The different 
tints denote regions with opposite sign. 



Although the matrix elements of the Coulomb potential 
are formally spin-dependent due to the spin-dependence 
of the MLWFs, we find that this dependence is negligible 
in practice, i.e., ~ ~ W^. 

In Fig. inja) we present a convergence study for the av- 
erage on-site diagonal matrix elements of the bare and 
the screened Coulomb interaction (U) between the 3d 
orbitals as a function of the number of k points for the 
NM and the FM states of bcc Fe. Fig.^b) shows the 
same matrix elements as a function of the number of 
bands used in the construction of the MLWFs. We ob- 
serve a fast k-point convergence but a relatively slow con- 
vergence with respect to the number of bands. In fact, 
the screened U is not completely converged even with 60 
bands. The large difference between the NM and FM 
states will be discussed below. 

The increase of the matrix elements of the Coulomb 
potential with the number of bands can be explained by 
the localization of the Wannier functions. In Fig.|^c) we 
show the spread f2 for the 3d orbitals as a function of the 
number of bands. It can clearly be seen that f2 decreases 
if the number of bands increases, indicating that the 3d 
orbitals become more localized, which in turn gives rise 
to a larger U . For the rest of the 3d transition-metal 
series the behavior of the Coulomb matrix elements with 
respect to the number of bands is very similar to bcc 
Fe. In the rest of this section we use 6 bands and an 
8x8x8 k-point mesh in order to compare our results 
with previously published data that adopted the same 
parameter settings. 

In TableUwe present the on-site matrix elements of the 
screened Coulomb potential for the NM and FM states 
of bcc Fe. TablesHIl and IIIII contain the values for the FM 
state of fee Co and Ni. For all three systems the aver- 
age values for the diagonal (U) and off-diagonal ([/', J) 
matrix elements are given in Table lTVl We note that the 
splitting of the diagonal matrix elements by the crystal- 
field effect is quite pronounced. For bcc Fe the e g -like 
diagonal elements are larger than the t2g-hke ones in the 
FM state, while in the NM state it is just the opposite. 
This means that the strongest interaction takes place be- 
tween the electrons in the e g -like (tag-like) orbitals for 
the FM (NM) state, because these are more localized. In 
the case of fee Ni the situation is very similar. However, 
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FIG. 3: (Color online) Convergence of the average on-site 
matrix elements of the bare and screened Coulomb interaction 
between the 3d orbitals as a function of (a) the number of k 
points and (b) the number of bands used in the construction 
of the MLWF basis for the NM and FM states of bcc Fe. (c) 
The same as (b) for the average spread of the 3d orbitals. 



in fee Co the splitting of the diagonal Coulomb matrix el- 
ements by the crystal- field effect is different: In the NM 
state (results not shown) all diagonal elements assume 
similar values, while in the FM state (see TableHI|) the 
t2 g -like diagonal elements are larger than the e g -like ones. 

Figure^] shows the average bare and screened on-site 
direct (U) Coulomb matrix elements for the series of 3d 
transition metals in the NM state. Results for the FM 
state of Fe, Co, and Ni are also included. Note that 
among all considered systems only the first three ele- 
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FIG. 4: (Color online) Average bare and screened on-site di- 
rect (U) Coulomb matrix elements between the 3d orbitals 
for the series of 3d transition metals. For comparison results 
from Ref.[13 (filled spheres) are given. We include 6 bands in 
the construction of the MLWFs. 



ments are not magnetic, while Cr and Mn order antifcrro- 
magnetically and Fe, Co, and Ni are ferromagnetic. The 
bare U for the NM state increases linearly from 14 eV for 
Sc to 25 eV for Ni. This stems from the fact that, as one 
moves from the left to the right within one row of the pe- 
riodic table, the nuclear charge increases and causes the 
3d wave functions to contract. Hence the localization 
of the 3d electrons increases, giving rise to the observed 
trend for U. However, this trend is not observed for 
the screened Coulomb interaction, where the calculated 
values lie between 0.8 eV and 1.5 eV. As already seen 
in TableU the matrix elements of the screened Coulomb 
potential depend on the magnetic state. Figured] indi- 
cates that the U values for the NM and FM states of Fe, 
Co, and Ni are also rather different, and this difference 
increases with the exchange splitting in the Ni-Co-Fe se- 
quence. This observation can be qualitatively explained 
by the density of states (DOS) around the Fermi level 
presented in Fig.[SJ As the screened Coulomb interaction 
depends on the polarizability, the number of occupied 
and unoccupied states around the Fermi level plays an 
important role in determining its strength. Bcc Fe in the 
NM state has the largest DOS around the Fermi energy 
and hence the smallest Coulomb matrix elements. How- 
ever, for FM Fe the majority and minority-spin peaks 
at the Fermi level are shifted to lower and higher ener- 
gies, respectively, due to the exchange field, leading to a 
lower DOS at the Fermi level. As a consequence, we ob- 
tain larger matrix elements. For the last three elements 
the calculated U, U', and J values increase linearly for 
the NM state (see Table HV]) . while they are almost con- 
stant for the FM state. A comparison of our results for 
the screened Coulomb interaction in the NM state with 
Ref.H3 is given in Fig.SJ The agreement for the U values 
is very good, and we find an equally good agreement for 
the J values, which are not displayed in the figure. 
To investigate the effect of static correlation on the 
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TABLE I: Screened on-site direct (U m n = Wmm,,nn) and exchange (J m „ = VFmS, n m) Coulomb matrix elements between the 3d 
orbitals for FM bcc Fe within LSDA. In parentheses we show results for the NM state. We include 6 bands in the construction 



of the MLWFs. The indices 1 and 2 (3, 4, and 5) correspond to the e g -like (£2 g - 


-like) Wannier orbitals. All enerj 


;ies are in eV. 


U mn 


1 


2 


3 


4 


5 


l 


1.63 (0.62) 


0.35 (0.05) 


0.64 (0.20) 


0.64 (0.20) 


0.31 (0.07) 


2 


0.35 (0.05) 


1.63 (0.62) 


0.42 (0.12) 


0.42 (0.12) 


0.75 (0.25) 


3 


0.64 (0.20) 


0.42 (0.12) 


1.33 (0.96) 


0.38 (0.17) 


0.38 (0.17) 


4 


0.64 (0.20) 


0.42 (0.12) 


0.38 (0.17) 


1.33 (0.96) 


0.38 (0.17) 


5 


0.31 (0.07) 


0.75 (0.25) 


0.38 (0.17) 


0.38 (0.17) 


1.33 (0.96) 


J mn 


1 


2 


3 


4 


5 


l 




0.64 (0.28) 


0.41 (0.30) 


0.41 (0.30) 


0.56 (0.37) 


4 


0.64 (0.28) 




0.51 (0.35) 


0.51 (0.35 


0.35 (0.27) 


3 


0.41 (0.30) 


0.51 (0.35) 




0.47 (0.41) 


0.47 (0.41) 


2 


0.41 (0.30) 


0.51 (0.35) 


0.47 (0.41) 




0.47 (0.41) 


5 


0.56 (0.37) 


0.35 (0.27) 


0.47 (0.41) 


0.47 (0.41) 





TABLE II: The same as Table|I] for the screened Coulomb 
potential for the FM state of fee Co. 



Umn 


1 


2 


3 


4 


5 


l 


1.20 


0.22 


0.54 


0.51 


0.24 


2 


0.22 


1.20 


0.36 


0.32 


0.59 


3 


0.54 


0.36 


1.45 


0.38 


0.39 


4 


0.51 


0.32 


0.38 


1.45 


0.36 


5 


0.24 


0.59 


0.39 


0.36 


1.45 


J mn 


1 


2 


3 


4 


5 


l 




0.50 


0.42 


0.41 


0.56 


2 


0.50 




0.53 


0.50 


0.36 


3 


0.42 


0.53 




0.54 


0.53 


4 


0.41 


0.50 


0.54 




0.55 


5 


0.56 


0.36 


0.53 


0.55 





TABLE IV: Average screened on-site direct (diagonal U and 
off-diagonal U') and exchange (J) Coulomb matrix elements 
between the 3d orbitals for FM 3d transition metals within 
LSDA. In parentheses we show results for the NM states. For 
comparison the LSDA+[7 results (U = 1.9 eV, J = 1.2 eV) for 
fee Ni are presented. We include 6 bands in the construction 
of the MLWFs. All energies are in eV. 





V 


& 


J 


bcc Fe 


1.45 (0.82) 


0.47 (0.15) 


0.48 (0.34) 


fee Co 


1.35 (1.04) 


0.39 (0.23) 


0.49 (0.40) 


fee Ni 


1.41 (1.28) 


0.41 (0.34) 


0.50 (0.46) 


fee Ni* 


1.49 


0.44 


0.51 



LSDA+C/ 



TABLE III: The same as Tablc[ll]for fcc Ni. 















Umn 


1 


2 


3 


4 


5 


l 


1.53 


0.28 


0.57 


0.58 


0.26 


2 


0.28 


1.53 


0.35 


0.37 


0.66 


3 


0.57 


0.35 


1.33 


0.34 


0.33 


4 


0.58 


0.37 


0.34 


1.33 


0.33 


5 


0.26 


0.66 


0.33 


0.33 


1.33 


J mn 


1 


2 


3 


4 


5 


l 




0.63 


0.43 


0.43 


0.58 


2 


0.63 




0.52 


0.53 


0.37 


3 


0.43 


0.52 




0.49 


0.49 


4 


0.43 


0.53 


0.49 




0.49 


5 


0.58 


0.37 


0.49 


0.49 





servation can be explained by the scenario given above. 
Within the LSDA+£7 scheme the exchange splitting of 
the Ni 3d states increases. This in turn gives rise to 
a larger magnetic moment (see TablelVT) and a reduced 
DOS around the Fermi level. As a consequence, the 
Coulomb matrix elements increase. We expect a similar 
behavior for bcc Fe and fcc Co if the LSDA+C/ scheme 
is employed. 

Finally, we discuss the values of the matrix elements 
of the screened Coulomb potential W used in the calcu- 
lation of R to make a connection with next section. The 
magnetic response function can be schematically written 
as 

K 

R = , (33) 

1 - WK ' V ' 



matrix elements of the screened Coulomb interaction we 
present LSDA+C/ results for U, U', and J for the FM 
state of Ni in Table llVl The average matrix elements 
slightly increase relative to the LSDA. Again this ob- 



where the screened Coulomb potential W in the denom- 
inator is responsible for the formation of collective spin- 
wave excitations. As shown above, the matrix elements 
of W depend on the number of bands included in the con- 
struction of the MLWFs. The matrix elements of K are 
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FIG. 5: (Color online) Spin-resolved density of states for bcc 
Fe, fee Co, and fee Ni for the NM and FM states. The vertical 
dashed lines denote the Fermi level. Note that positive val- 
ues of DOS refer to the majority-spin electrons and negative 
values to the minority-spin electrons. 



considerably less sensitive to the number of bands; usu- 
ally 10 bands are sufficient for 3d ferromagnets. For this 
reason, in the calculation of R one might choose a partic- 
ular W that satisfies the exact condition lim q _ ) .o w(q) = 
(Goldstone mode). For instance, in fee Ni within LSDA 
one should include about 100 bands in order to fulfill the 
Goldstone theorem. Alternatively, one can calculate W 
for a given number of bands and then scale it by a factor 
?7, i.e., W — > r/W, to obtain the Goldstone mode correctly. 
This second approach is computationally less demanding 
and is used in the present work. For 3d ferromagnets 
we calculate W by including only 6 bands. In the case 
of fee Ni the scaling factor is r\ w 1.5 within LSDA and 
7/ w 1.8 within LSDA+f/. This means that the Coulomb 
matrix elements presented in Table llLTl should be multi- 
plied by 1.5 in the calculation of the magnetic response 
function R within LSDA for fee Ni. The violation of the 
Goldstone theorem within the present formalism stems 
from the approximations made in the calculation of the 
kernel K and the screened Coulomb potential W . In 
a fully self-consistent linear-response calculation without 
additional approximations the Goldstone theorem should 
be fulfilled. However, this is hardly feasible in practice. 
When we manually reduce the exchange splitting of Ni 
by one-half within LSDA to simulate the renormalized 
Green function, the scaling factor is reduced to 77 ~ 1.1. 
Note that the LSDA overestimates the exchange splitting 
of Ni by a factor of 2 compared to experiments / ' 71 In 
systems like bcc Fe for which the LSDA already provides 
a reasonable description of the electronic band structure 
compared to the experiment a 71 i 73 i 74 the scaling factor ij 
is close to 1. 



TABLE V: Calculated magnetic moment m, exchange split- 
ting at X2 and L3 as well as spin-wave stiffness constant D 
in fee Ni within LSDA and LSDA+{7 with U = 1.9 eV and 
J = 1.2 eV. For comparison experimental values are given. 





m(p B ) 


X 2 (eV) 


L 3 (eV) 


D(meVA 2 ) 


LSDA 


0.61 


0.61 


0.57 


740 


LSDA+t/ 


0.65 


0.55 


0.62 


540 


Expt. 


0.60 a 


Q.20 b 


0.30 c 


550 d 



a Ref. Q 
^Ref. 7(1 
c Ref. 7i| 
d Ref. [H 



IV. MAGNETIC EXCITATIONS IN fee Ni 

This section deals with magnetic excitations in fee Ni. 
Among the 3d ferromagnets Ni is known for particu- 
larly large discrepancies between the results from DFT 
calculations and experiments: The width of the occu- 
pied 3d bands in the LSDA is about 30% larger than 
that found in photoemission experiments, whereas the sp 
band width agrees within 10%^ 7 -^ Similarly, the LSDA 
yields a much smaller DOS [1.9 states/ (eV atom)] at 
the Fermi level compared to low-temperature specific- 
heat data [3.0 states/(eV atom)], indicating a quasipar- 
ticle mass enhancement^ Even larger discrepancies are 
obtained for the exchange splitting. Photoemission ex- 
periments give a small and highly anisotropic exchange 
splitting, 0.3 eV at the L3 point and 0.2 eV at the X2 
point i 70 ' 71 In contrast, the LSDA yields a rather large 
(0.6 eV) and almost isotropic splitting. 67 However, the 
calculated magnetic moment turns out to be in good 
agreement with the experimental valued 

Our calculated magnetic moments and exchange split- 
tings within the LSDA and LSDA+C7 are presented in 
TablelVl The corresponding band structures are given in 
Fig. [5] Within the LSDA the band structure, exchange 
splitting, and magnetic moment are in good accordance 
with literature valuesi 67 ' 68 ' 72 The inclusion of explicit 
static correlation in the form of a Hubbard U within 
LSDA+J7 slightly changes the electronic structure of Ni. 
In this case the exchange splitting is more anisotropic 
compared to the LSDA, i.e., it increases at the L3 point 
and decreases at the X 2 point (see Table lV|) . The average 
exchange splitting increases within the LSDA+i7 scheme 
however, and depending on the values of the Hubbard 
parameters U and J the LSDA+C/ hence gives rise to a 
larger magnetic moment. 

The magnetic response function R [see Eq. ([3"3")l ] con- 
tains all relevant information about the dynamics of the 
spin system. The poles of K in the numerator corre- 
spond to the energies of single-particle spin-flip Stoncr 
excitations, while the zeros of the denominator (1 — WK) 
describe collective spin-wave excitations. In the preced- 
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FIG. 6: (Color online) LSDA and LSDA+(7 (U = 1.9 eV and 
J = 1.2 eV) band structures of FM fee Ni for (a) majority 
and (b) minority spins. 
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FIG. 7: (Color online) The Kohn-Sham magnetic response 
function for fee Ni for selected wave vectors along the T — X 
direction in the Brillouin zone. Since the Kohn-Sham system 
is noninteracting, the spectrum exhibits only single-particle 
Stoner excitations but no collective magnon modes at lower 
energies. In each panel the peak amplitudes are scaled to the 
same height. 



ing section we showed that the convergence of the matrix 
elements of the screened Coulomb potential W with re- 
spect to the number of k points (number of bands) is 
fast (slow), whereas the situation is opposite for the ker- 
nel K. For W we take the values from Sec. lIIIl scaled 
by an appropriate factor 77, while we construct K using 
5 MLWFs, 15 occupied and unoccupied bands per spin 
channel, and a very dense 40 x 40 x 40 k-point sampling. 
The Brillouin-zone summations in K are performed with 
the tetrahedron method^ The calculation of R is car- 
ried out for a fixed q as a function of energy uj up to 1.5 
eV. 

The remainder of this section is divided into two parts. 
In the first part the single-particle Stoner excitations are 
presented. The second part deals with the collective spin- 
wave excitations. 



A. Single-particle Stoner excitations 

Stoner excitations are electron transitions between 
bands of opposite spin. When an electron is excited from 
an occupied majority-spin state at k to an unoccupied 
minority-spin state at k+q, it produces an electron-hole 
pair with triplet spin configuration that reduces the mag- 
netization by unity. Therefore, these excitations are as- 
sociated with longitudinal fluctuations of the magneti- 
zation and play an important role in determining the 



high-temperature properties of magnetic materials^ To 
simplify the discussion, let us consider a free electron gas 
with parabolic dispersion and a rigid exchange splitting 
AE ex . The corresponding single-particle excitation ener- 
gies are given by w(q) = (±2k-q+q 2 ) + AE CX . For q — > 
the necessary energy for such excitations equals the ex- 
change splitting AE CX , but for finite q these excitations 
form a continuous spectrum, which is called the Stoner 
continuum and is determined by all possible values of k. 
The upper and lower bounds of the Stoner continuum 
depend on the electronic band structure of the material. 
For example, in strong ferromagnets the smallest possi- 
ble excitation energy is given by the Stoner gap AS. On 
the other hand, weak ferromagnets exhibit a vanishing 
gap AS — 0, and thus single- particle excitations do not 
require a finite energy for particular q values. The lower 
bound is of special interest, because when the collective 
spin-wave excitations enter the Stoner continuum, they 
start to decay into single-particle excitations, which re- 
duces their lifetime drastically. We note that the picture 
of Stoner excitations in real materials is different from 
the free electron gas with a single band^S 

Since the Stoner excitations are single-particle spin-flip 
processes, they can be studied qualitatively at the Kohn- 
Sham level. They result from the term Im K, as discussed 
before. In Fig. [7] we present the imaginary part of the 
Kohn-Sham magnetic response function for fee Ni for se- 
lected wave vectors along the T — X direction. Note that 
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all peak amplitudes are scaled to the same height. The 
energetic position and shape of the Stoner spectrum for 
q = provide a measure for the mean exchange splitting 
and its variation across the Brillouin zone. The triple- 
peak structure reflects the intraband and interband tran- 
sitions. The position of the first peak gives the mean 
exchange splitting: 0.61 eV within LSDA and 0.69 eV 
within LSDA+[7. The broadening of the peaks reflects 
the k-dependence of the exchange splitting. For rigidly 
split bands one would get S peaks. As the wave vector 
increases, individual peaks become broader and are no 
longer distinguishable. Since the Kohn-Sham system is 
noninteracting, the spectrum does not exhibit collective 
magnon modes at low energies. 

Experimentally, Stoner excitations in id ferromag- 
nets are studied with spin-polarized electron-energy-loss 
spectroscopy (SPEELS), a technique that not only mea- 
sures the high-energy Stoner excitations but also low- 
energy collective spin-wave modes up to the Brillouin- 
zone boundary. Using SPEELS Kirschner et alS^ stud- 
ied the Stoner excitation spectrum of Ni at q ~ 0. The 
authors found that the spectrum had a broad energy dis- 
tribution of 0.3 eV (full width at half maximum) and 
was centered around 0.3 eV, which is consistent with the 
average exchange splitting determined by photoemission 
experiments. Additionally, the width of the distribu- 
tion provided evidence for the pronounced k-dependence 
of the exchange splitting. The comparison of our cal- 
culated spectra with the experimental data shows large 
discrepancies, however. These can be attributed to the 
overestimated exchange splitting and its incorrect nearly 
isotropic behavior over the Brillouin zone in the LSDA. 
As pointed out by Oles and Stollhoff 78 and by Liebsch, 79 
the large exchange splitting in the LSDA is due to the 
neglect of strong correlation effects within the 3d states 
and anisotropic exchange. Consequently, recent studies 
by Katsnelson and Lichtensteir* 8 -^ and by Grechnev et 
al£^ based on dynamical mean-field theory, in which the 
exchange and correlation effects are taken into account 
properly, gave much improved results for the exchange 
splitting and quasiparticle band structure of Ni. 



B. Collective spin- wave excitations 

In the preceding section we discussed the non- 
interacting magnetic response function K, which has sin- 
gularities that correspond to single-particle Stoner exci- 
tations, usually with a broad energy distribution. In ad- 
dition, there may be other singularities of R for a fixed q. 
For some w(q) outside the Stoner continuum the denom- 
inator (1 — WK) of the magnetic response function might 
vanish, indicating collective spin-wave excitations. These 
correspond to transverse fluctuations of the direction of 
the magnetization and can be interpreted as a coherent 
superposition of an electron and a hole with opposite 
spins coupled via an attractive screened Coulomb inter- 
action W, thereby forming a bound state with energy 



a;(q), spin 1, and momentum q^ As more than one par- 
ticle is involved in the excitation process, the formation 
of spin waves cannot be described within a single-particle 
picture. 

Before discussing our numerical results we briefly re- 
view the status of experimental and theoretical studies 
of magnetic excitations in fee Ni. Starting from the 
mid 1960s the spin dynamics of ferromagnetic 3<i tran- 
sition metals and their alloys were intensively investi- 
gated by inelastic neutron scattering. The first exper- 
iments for Ni were performed at high temperatures and 
in the low-energy region. Later the measurements were 
extended to higher energies. The results of these early 
neutron-scattering experiments and a comparison with 
realistic band structure calculations can be found in the 
review article by Lowde and Windsor^ However, the 
signal intensity and energy resolution in these early ex- 
periments were not favorable for a quantitative deter- 
mination of the spin-wave excitations in Ni. More pre- 
cise measurements were reported by Mook and collabo- 
rators in the 1980s with emphasis on high energies up 
to 240 meVi^— The authors measured the spin-wave 
dispersion of Ni at several temperatures starting from 
T = 4.2 K up to T « 2T C , where T c = 631 K is the 
Curie temperature, and found that the spin- wave disper- 
sion was isotropic in q over the entire temperature range 
studied. The obtained spin-wave stiffness constant was 
D = 550 meV A 2 at T = 4.2 K and D = 505 meV A 2 
at T = 295 K. The spin-wave intensity was found to de- 
crease faster with increasing energy in the [1 1 1] direction 
than along other symmetry directions. In all directions 
the spin waves eventually disappeared at some wave vec- 
tor close to the zone boundary inside the first Brillouin 
zone, which was attributed to a decay into Stoner exci- 
tations. Additionally, the authors found evidence for a 
second branch in the spin- wave dispersion, i.e., an optical 
mode that crosses the main [1 0] acoustic branch around 
125 meV. However, inelastic neutron scattering does not 
allow to explore the entire Brillouin zone in 3d ferromag- 
nets, in contrast to SPEELS. Using SPEELS Abraham 
and Hopster 87 attempted to study short-wave-length (or 
large-wave-vector) spin excitations in Ni. However, they 
only detected Stoner excitations, although the data re- 
ported reaches down to 100 meV and the resolution of 
the instrument (17 meV) should have been sufficient to 
observe collective excitations. A qualitative explanation 
for the absence of spin-wave peaks in the SPEELS spec- 
trum was given by Hong and Mills ^ who showed that 
the spin waves can only be observed in SPEELS if the 
exchange splitting of the 3d bands is large compared to 
the spin-wave excitation energies. Indeed, spin waves up 
to the Brillouin-zone boundary are observed in Fe and 
Co, which have substantially larger exchange splittings 
than NhS 8 .^ 

On the theoretical side, the first attempt to use real- 
istic energy bands for Ni in the calculations of a gen- 
eralized susceptibility was undertaken by Lowde and 
Windsor. 8 - The authors calculated the magnetic suscep- 
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tibility within the RPA for rigidly spin-split bands, but 
the agreement with the available experimental data was 
not good. Cooke et a/.— ~— showed that it is necessary 
to take the k-dependence of the exchange splitting into 
account for a quantitative comparison between theory 
and experiment. They used a tight-binding description 
of the electronic energy bands, and the ferromagnetism 
was driven by an empirical on-site Coulomb interaction 
between the 3d electrons with two adjustable parame- 
ters. These parameters were chosen in such a way that 
the calculations reproduce the experimentally observed 
magnetic moment as well as the correct t2 g and e g char- 
acter of the moment as measured in neutron magnctic- 
form-factor experiments. The calculations of Cooke et 
at not only yielded the correct spin-wave dispersion re- 
lation, including the appearance of the optical branch in 
the [1 0] direction, but the damping of the spin waves 
in the presence of the Stoner modes was also correctly 
described. A similar approach was used by Hong and 
Mills 92 with an empirical Coulomb interaction that was 
form-invariant under spin rotations. However, they failed 
to find the optical mode in the spin-wave dispersion of 
Ni. 

A much more accurate description of spin waves in 
ferromagnetic 3d transition metals was reported by 
Savrasov^i and by Karlsson and Aryasetiawan* 3 ^, In both 
works spin-polarized DFT was used for the ground- 
state calculations. For the transverse spin susceptibility 
Savrasov employed TDDFT, while Karlsson and Aryase- 
tiawan adopted MBPT. Similar to the findings of Cooke 
et at, Savrasov obtained two branches in the spin- wave 
dispersion of Ni along the [1 0] direction with an opti- 
cal mode at high energies. Karlsson and Aryasetiawan 
confirmed these results and investigated the role of the 
one-particle band structure on the spin-wave dispersion. 
They found a very good agreement between theory and 
experiment if the LSDA exchange splitting was manually 
reduced by one-half. Additionally, the calculations gave 
evidence for an optical branch along the [1 1 1] direction. 

The optical branch in the spin- wave dispersion of Ni is 
evidently a very subtle issue. So far, there has been no 
general consensus in theoretical treatments concerning 
the sensitivity of the results to the details of the elec- 
tronic structure and to the method used. In the fol- 
lowing we hence focus on the energy region where the 
double-peak structure is observed in the spin- wave spec- 
trum. To illuminate the effect of the electronic structure 
on the magnetic excitation spectrum of Ni we employ 
three different methods: LSDA, LSDA+C7, and LSDA 
with a reduced exchange splitting by one-half. The cal- 
culated spin-wave dispersion along the high-symmetry 
lines L — r — X is displayed in Fig. [5] For comparison 
the experimental dispersion is also shown. The LSDA 
and LSDA+J7 yield qualitatively similar results, with an 
optical mode not only in the [1 0] but also in the [1 1 1] 
direction. The acoustic branch is well described within 
LSDA and LSDA+J7, but the optical branch is too high 
in energy. This discrepancy between theory and exper- 
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FIG. 8: Spin-wave dispersion for fee Ni along high-symmetry 
lines (L — r — X) in the Brillouin zone within (a) LSDA, (b) 
LSDA+C/, and (c) LSDA with a reduced exchange splitting. 
The experimental dispersion is taken from Refs.Hg and l86l. 



iment can be traced back to the overestimation of the 
exchange splitting in the LSDAi22 Indeed, when we re- 
duce the exchange splitting by one-half as in Ref.[33l we 
obtain reasonable agreement with the experiments. The 
corresponding dispersion, shown in Fig.[Stc) , is similar to 
that of Karlsson and Aryasetiawan.— 

In Fig. [9] we show the imaginary part of the magnetic 
response function R for selected wave vectors along the 
[1 0] and [1 1 1] directions. The peak amplitudes are 
again scaled to the same height for presentational pur- 
poses. The obtained spin- wave spectra along [1 0] are 
in very good agreement with previous calculations! 31 ! 33 
In particular, a double-peak structure starts to develop 
from q = (0.15, 0, 0). For larger wave vectors the two 
peaks overlap, resulting in a rather broad single feature, 
which can be decomposed into two Lorentzian peaks as 
shown in Fig. [TU] for q = (0.25, 0, 0). As the wave vec- 
tor increases, the intensity of the lower peak decreases 
in agreement with experiments.— In the [1 1 1] direction 
the double peak is not so clear within the LSDA. Never- 
theless, the calculated structure can still be decomposed 
into two Lorentzians (results not shown), but the q in- 
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FIG. 9: (Color online) Spin-wave excitation spectra for fee Ni 
within LSDA and LSDA+C7 for selected wave vectors along 
(a) the [1 0] and (b) the [1 1 1] directions in the Brillouin 
zone. The double-peak structure for certain wave vectors is 
clearly visible. All peak amplitudes are scaled to the same 
height. 



terval where the double-peak structure appears is small 
compared to the [1 0] direction. 

Our calculations show that the exchange splitting has a 
strong influence on the emergence of a double-peak struc- 
ture. To demonstrate this we consider two q points in the 
[1 0] direction where such features appear and calculate 
the LSDA spin-wave spectra by reducing the exchange 
splitting gradually by 0.1 eV in each step, up to one-half 
of its original value. The results are presented in Fig.lTTl 
Even a very small reduction in the exchange splitting by 
0.1 eV strongly suppresses the double-peak feature. As 
the exchange splitting reduces further, the width of the 
peaks becomes narrower, indicating an increased spin- 
wave life time. When the exchange splitting is reduced 
by one-half and corresponds to the experimental value, 
the optical branch completely disappears in the [1 1 1] 
direction, while it is shifted to larger wave vectors in 
the [1 0] direction. The situation is very similar within 
LSDA+J7 with a reduced exchange splitting. Karlsson 
and Aryasetiawan found a weak double-peak structure 
with a reduced exchange splitting for q = 0.1875(1, 1, 1) 
when they used a small Gaussian broadening parameter 
in the Brillouin-zone integration^ but this double-peak 
structure becomes smeared out if a larger broadening is 
used. Based on the symmetry properties of fee Ni, the 




ra(eV) 

FIG. 10: (Color online) Decomposition of the single feature 
for q = (0.25, 0, 0) into two Lorentzian peaks. In the inset we 
compare with the original curve. 



spin- wave dispersion should be isotropic, i.e., one expects 
an optical branch also in the [1 1 1] direction, but the reso- 
lution of the presently available experimental data is not 
sufficient to observe iti^ Of course, a reduction of the 
exchange splitting in the calculations does not solve all 
problems for Ni. In particular, too large band widths and 
the incorrect isotropic exchange splitting remain impor- 
tant issues whose effect on the spin- wave spectra requires 
further investigations. 

As pointed out by Cooke et oZ.j=« the double-peak 
structure in the magnetic excitation spectra of 3d ferro- 
magnets stems from the k-dependent exchange splitting 
and interband transitions. A detailed analysis was given 
by Karlsson and Aryasetiawan^ who showed that the 
double-peak structure in fee Ni is implicitly contained 
in the kernel K, i.e., it is a band-structure effect arising 
from the fact that WIm K possesses additional structure 
below the Stoner peak, which in turn gives rise to a weak 
dip structure in (1 — WKeK) via the Hilbert transfor- 
mation. Such a weak dip structure appears as a second 
peak in the spin- wave excitation spectrum for particular 
wave vectors. This observation is indirectly confirmed 
by our analysis. The application of a Hubbard U in the 
LSDA+f/ calculation shifts majority-spin and minority- 
spin states more or less rigidly, and the optical branch 
appears in the LSDA and LSDA+[/ around the same 
wave vector q along [10 0], while a reduction of the ex- 
change splitting leads to an appearance of the optical 
branch at higher q values. 

The magnetic response function R allows to extract 
information about the spin- wave life times, which are in- 
versely proportional to the widths of the spin- wave peaks. 
If Ni were a strong ferromagnet, one would get well de- 
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FIG. 11: (Color online) Same as Fig.[9]for q = (0.15, 0, 0) and 
q = (0.2, 0, 0) for different exchange splitting within LSDA. 
The spectra for the original LSDA exchange splitting are in- 
cluded for comparison. 

fined S peaks up to about 0.3 eV in the spin-wave spec- 
tra, corresponding to the energy difference between the 
highest occupied majority-spin 3d band and the Fermi 
level, i.e., the Stoner gap. However, Ni is not a truly 
strong ferromagnet due to the sp-d hybridized majority 
states around the Fermi energy, and Stoner excitations 
thus occur essentially at all energies. This means that 
spin waves decay into Stoner excitations for all non-zero 
wave vectors, which is reflected by a finite width of the 
spin- wave peaks as shown in Fig. [HI The life time of the 
spin waves depends on the details of the coupling be- 
tween these excitations and on the density of states of 
the Stoner excitations. It should be noted that in itiner- 
ant ferromagnets the life time of the spin waves becomes 
infinite in the limit of a vanishing wave vector, provided 
that spin-orbit coupling is ignored. As seen in Fig.[9j for 
small wave vectors the spin-wave peaks are indeed nar- 
row, and the damping is hence weak, i.e., the life times 
are long. As the wave vector increases, the spin- wave 
dispersion enters into the region with a high density of 
states of Stoner excitations, and the decay mechanism 
becomes more efficient. In both crystallographic direc- 
tions considered here the maximum damping occurs in 
the region where the double-peak structure appears. In 
contrast to the findings of Cooke et al.^ the spin- wave 
peaks associated with the optical branch are narrower 
than the acoustic ones in our work, and the peak widths 
stay almost constant throughout the Brillouin zone, while 
the peak widths in the acoustic branch increase with the 
wave vector. This might in fact explain why the acoustic 
branch disappears in the middle of the Brillouin zone in 
the neutron-scattering experiments^ 

Finally, we focus on the spin- wave stiffness constant D. 
In ferromagnets the spin waves show a quadratic disper- 



sion law Lu(q) — Dq 2 for small wave vectors. The values 
for Ni obtained with different methods are listed in Ta- 
ble[V] Our LSDA estimate of 740 meV A is substantially 
larger than the experimental value 550 meV A . With the 
reduced exchange splitting D increases even further to 
870meVA , whereas LSDA+C/ provides a much better 
estimate of 540meVA , which reflects the importance 
of static correlation effects in Ni. We note that our 
LSDA estimate for the spin-wave stiffness of fee Ni is in 
good agreement with calculations based on constrained 
DFT in the adiabatic approximation. Using the frozen- 
magnon technique Rosengaard and Johansson 18 found 
D = 739meVA , which is very similar to the values ob- 
tained by Schilfgaarde and Antropov^ (740 meV A ) and 

by Pajda et al££ (756 ± 29meVA ), who employed real- 
space methods. As the adiabatic approximation becomes 
exact in the limit of long wave lengths (q — > 0),— this 
can be compared with values obtained from more rigor- 
ous approaches based on the dynamical transverse spin 
susceptibility or magnetic response function. However, 
within constrained DFT the second branch in the spin- 
wave dispersion of Ni cannot be described. 



V. CONCLUSIONS AND OUTLOOK 

In conclusion, we have developed a computational 
method to study excitation spectra of magnetic materi- 
als from first principles. The method is based on many- 
body perturbation theory. The main quantity of inter- 
est is the transverse magnetic response function, which 
treats both collective spin-wave excitations and single- 
particle spin-flip Stoner excitations on an equal footing. 
In order to describe the former we include appropriate 
vertex corrections in the form of a multiple-scattering 
T-matrix, which describes the coupling of electrons and 
holes with different spins. To reduce the numerical cost 
for the calculation of the four-point T-matrix we exploit 
a transformation to maximally localized Wannier func- 
tions that takes advantage of the short spatial range of 
electronic correlation in the partially filled d or / orbitals 
of magnetic materials. Our implementation is based on 
the FLAPW method. 

The developed scheme was employed to calculate the 
matrix elements of the Coulomb potential for the series 
of 3c? transition metals in the MLWF basis. Special at- 
tention was given to the ferromagnets Fe, Co and Ni. We 
showed that the matrix elements of the screened Coulomb 
potential are rather different for the NM and FM states 
and that the difference increases with the exchange split- 
ting in the Ni-Co-Fe sequence, which can be accounted 
for on the basis of the total density of states around the 
Fermi level for the corresponding systems. 

The magnetic excitations in fee Ni were studied in de- 
tail based on the LSDA and LSDA+C7 methods. Both 
schemes give qualitatively similar results for the spin- 
wave spectra and dispersion. However, correlation effects 
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seem to be important for the spin- wave stiffness constant, 
which is overestimated within LSDA. Our calculations 
indicate the existence of an optical branch in the spin- 
wave dispersion of Ni along the [1 1 1] in addition to that 
along the [1 0] direction in the Brillouin zone. Although 
the acoustic branch is well described within LSDA and 
LSDA+J7, the optical branch appears to be too high in 
energy. This discrepancy between theory and experiment 
can be attributed to the overestimation of the exchange 
splitting or, in other words, to the use of the Kohn-Sham 
Green function in the calculation of the kernel K instead 
of the renormalized one. 

In the LSDA Kohn-Sham Green function the long and 
short-range correlation effects arc not taken into account 
properly. The former can be treated within the GW 
approximation, while the latter require the summation 
of spin-dependent T-matrix contributions. In the future 
we plan to incorporate electron-electron (hole-hole) and 



electron-magnon scattering processes into the electronic 
self-energy by means of the T-matrix formalism, which 
improves the theoretical description of the quasiparticle 
band structure. In particular, it is expected to yield the 
correct exchange splitting in magnetic materials. 
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